#coding=utf-8

'''
任选思路进行积分： 单位球体积，不限于：

- 投点法
- 均匀抽样
- 重要抽样
'''

import numpy as np
import matplotlib.pyplot as plt

# 投点法
nsample = 100000

# 在 x = [-1,1], y=[-1,1], z=[-1,1] 这样一个 2x2x2 的立方体里随机投点
x = 2*np.random.random(nsample) -1
y = 2*np.random.random(nsample) -1
z = 2*np.random.random(nsample) -1
count=0 # 统计有多少点在单位球内部
for i in range(nsample):
    if( x[i]*x[i] + y[i]*y[i] + z[i]*z[i] < 1 ):
        count += 1

mu = count / nsample; print("mu = ", mu)
sigma2 = ( mu - mu*mu )/nsample

print( " 5sigma result: volume ~ ", 8*np.array([mu - 5*np.sqrt(sigma2), mu + 5*np.sqrt(sigma2)]), " theory = ", 4/3*np.pi )

# 误差分析： 每次投点，投在球内的点数要么是１，要么是０，期望值是 mu, 方差是 mu - mu^2，方差可以用

# 蒙卡积分均匀取样: int^1_{-1} int^{ sqrt(1-x^2) }_{ - sqrt(1-x^2) } 2 sqrt(1-x^2 - y^2) dx dy
sumf = 0; sumf2 = 0; count = 0;
for i in range(nsample): # 在单位圆内均匀取点
    x = 2*np.random.random() -1
    y = 2*np.random.random() -1
    if x*x + y*y < 1:
        f = 2 * np.sqrt(1-x*x-y*y)
        sumf += f; sumf2 += f*f; count +=1
avef = sumf/count; avef2 = sumf2/count
mu = avef; sigma = np.sqrt((avef2 - avef*avef)/count)
print( " 5sigma result: volume ~ ", np.pi * np.array([ mu - 5*sigma, mu + 5*sigma ]) )

# 蒙卡积分重要抽样：省略